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1.  Introduction 


1.1  Objectives 

Results  from  previous  empirical  studies  indicate  the  high-frequency  Lg/P  ratio  is  one  of  the 
most  promising  discriminants  at  regional  distances.  However,  many  of  these  studies  are 
based  on  co-located  earthquakes  and  explosions  in  limited  geographic  regions,  and  their 
results  must  be  recalibrated  for  each  new  source  region.  Our  objective  is  to  develop  and 
apply. distance  corrections  for  up  to  10  stations  in  the  International  Data  Center  (IDC) 
primary  network  so  that  we  can  evaluate  the  transportability  of  the  Lg/P  discriminant  to 
uncalibrated  regions.  We  will  also  determine  the  sensitivity  of  the  Lg/P  discriminant  to  the 
accuracy  of  the  attenuation  models,  and  attempt  to  generalize  the  discriminant  by  geologic 
and  tectonic  regionalization. 

1.2  Background 

Several  early  studies  propose  the  use  of  Lg/P  ratios  as  a  regional  event  discriminant  for 
earthquakes  and  nuclear  explosions  [e.g.,  Willis  et  al,  1963;  Blandford,  1981;  Gupta  and 
Bumetti,  1981;  Nuttli,  1981;  Pomeroy  et  ai,  1982;  Murphy  and  Bennett,  1982;  Bennett  and 
Murphy,  1986;  Taylor  et  al,  1989].  These  studies  found  that  Lg/P  discriminants  provide 
some  separation  between  nuclear  explosions  and  earthquakes,  but  that  there  is  significant 
overlap  between  the  two  populations.  More  recent  studies  have  exploited  higher 
frequencies  to  improve  the  Lg/P  discriminant  and  have  extended  the  application  to  include 
industrial  chemical  explosions.  For  example,  high-frequency  (2-16  Hz)  Lg/P  ratios  have 
shown  to  be  successful  in  discriminating  between  earthquakes  and  mining  explosions  in 
northern  and  central  Europe  [e.g.,  Baumgardt  and  Young,  1990;  Dysart  and  Pulli,  1990; 
Baumgardt  et  al,  1992;  Wuster,  1993],  and  earthquakes  and  underground  nuclear 
explosions  in  Eurasia  [e.g.,  Bennett  et  al,  1989;  Chan  et  al,  1990;  Bennett  et  al,  1992]. 
The  practice  of  mining  can  also  induce  stress-release  events  such  as  rockbursts.  Often 
these  stress-release  events  are  much  larger  than  the  actual  mining  explosions.  Bennett  et  al. 
[1993]  show  that  stress-release  events  in  central  Europe  and  South  Africa  have  similar  Lg/ 
P  ratios  (>  1 .0)  to  earthquakes  across  broad  frequency  bands,  and  these  are  higher  than  the 
ratios  for  underground  nuclear  explosions  (<  1.0)  at  frequencies  above  2  Hz.  Most  of  these 
studies  use  co-located  explosions  and  earthquakes  to  minimize  propagation  effects.  A 
major  concern  expressed  in  the  recent  research  is  that  propagation  characteristics  may  have 
a  significant  affect  on  the  Lg/P  discriminant,  and  that  transportability  may  be  problematic 
[e.g,.  Lynnes  and  Baumstark,  1991;  Bennett  et  al.,  1992]. 

1.3  Current  Status 

This  report  describes  our  progress  in  the  first  year  of  this  two-year  effort.  We  worked  on 
compiling  data  sets  and  developing  software  for  generating  attenuation  models  and 
evaluating  the  Lg/P  discriminant.  We  modified  existing  software  to  measure  frequency- 
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dependent  amplitudes,  and  we  are  developing  new  software  to  invert  these  amplitudes  for 
source  and  attenuation  parameters. 

1.4  Outline  of  Report 

After  this  brief  introduction,  this  report  is  divided  into  three  main  sections.  Section  2 
describes  the  data  sets  needed  for  this  project  and  how  they  were  obtained.  Section  3 
describes  progress  on  our  attenuation  study  which  includes  investigating  alternative 
amplitude  measurements  and  inverting  them  for  source  and  attenuation  parameters.  We 
have  not  yet  started  the  discrimination  study.  Section  4  summarizes  our  progress  and 
describes  our  plans  for  the  second  year  of  this  project. 
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2.  Data  Sets 

This  section  describes  the  data  sets  that  we  are  collecting  to  support  our  attenuation  and 
discrimination  studies.  Our  primary  data  source  is  the  Group  of  Scientific  Experts  Third 
Technical  Test  (called  GSETT-3)  which  is  being  conducted  at  the  Prototype  International 
Data  Center  (PIDC)  in  Arlington,  Virginia  (for  an  overview  of  GSETT-3,  see  Kerr  [1993]). 
The  time  period  of  GSETT-3  is  from  January  1, 1995  to  the  present.  However,  we  limit  this 
report  to  include  data  prior  to  June  1,  1996.  Future  reports  will  include  more  recent  data. 

2.1  Station  Selection 

The  IDC  is  currently  receiving  data  from  46  primary  stations  and  87  auxiliary  stations.  As 
shown  in  Figure  1,  the  primary  stations  are  distributed  globally.  The  data  from  the  primary 
network  are  sent  continuously  to  the  IDC  where  they  are  processed  and  analyzed. 
Auxiliary  stations  expand  the  global  coverage  of  the  primary  network  (Figure  2),  but  these 
stations  are  used  only  when  data  are  requested  to  improve  the  location  solutions  for  events 
defined  by  the  primary  network.  We  select  stations  from  the  IDC  primary  and  auxiliary 
networks  to  include  in  our  study  based  on  the  ability  to  form  subnetworks  of  primary  and 
auxiliary  stations  that  are  in  close  proximity  (so  that  events  recorded  by  multiple  stations  at 
regional  distance  are  available  for  the  attenuation  study),  tectonic  regionalization,  and  the 
availability  of  regional  phases  in  the  IDC  Reviewed  Event  Bulletin,  or  REB. 

2.1.1  Subnetworks 

Most  primary  stations  have  other  primary  and  auxiliary  stations  in  close  proximity.  We 
give  preference  to  these  stations  because  events  recorded  at  multiple  stations  improve  the 
resolution  of  the  attenuation  models.  For  this  reason,  we  group  stations  into  subnetworks  if 
they  are  within  40°  of  each  other  (we  use  20°  as  the  largest  distance  possible  for  a  regional 
event).  Thus,  we  define  a  subnetwork  to  be  comprised  of  both  primary  and  auxiliary 
stations.  Some  IDC  primary  and  auxiliary  stations  may  not  be  included  because  they  did 
not  have  a  sufficient  number  of  regional  phases  (see  Section  2.1.3). 

2.1.2  Tectonic  Regionalization 

Eight  tectonic  regions  based  on  a  2°  by  2°  grid  are  defined  by  O.  Gudmundsson  of  the 
Australian  National  University  [personal  communication].  His  regionalization  is  shown  in 
Figure  3.  The  regions  are:  Young  Ocean  (<  25  Ma),  Intermediate  Ocean  (25-100  Ma),  Old 
Ocean  (>  100  Ma),  Tectonic,  Cenozoic  Continental,  Paleozoic  and  Mesozoic  (150-800 
Ma),  Proterozoic  (800-1700  Ma),  and  Archean  (>  1700  Ma).  The  Young  Ocean, 
Intermediate  Ocean,  and  Old  Ocean  categories  represent  the  age  of  the  ocean  floor.  The 
Tectonic  category  represents  young  and  tectonically-active  regions.  The  Cenozoic 
Continental,  Paleozoic  and  Mesozoic,  Proterozoic,  and  Archean  represent  young  to  old 
relatively  stable  continental  environments.  Figure  3  shows  the  distribution  of  IDC  primary 
stations  in  the  various  geologic  and  tectonic  environments.  Note  that  some  station 
locations  fall  on  the  boundary  of  two  adjacent  environments.  This  must  be  considered 
when  compiling  networks  that  represent  different  tectonic  regions.  For  the  purposes  of  this 
study,  we  combine  Cenozoic  Continental,  Paleozoic  and  Mesozoic,  Proterozoic,  and 
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Figure  1.  Location  of  IDC  primary  stations  participating  in  GSETT-3. 
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Figure  2.  Locations  of  IDC  auxiliary  stations  participating  in  GSETT-3. 


Figure  3.  Tectonic  regionalization  map  showing  locations  of  IDC  primary  stations.  This  tectonic  map  is  based 
2°  by  2°  grid  provided  by  O.  Gudmundsson  of  the  Australian  National  University  [personal  communication]. 


Al'chean  into  a  single  “stable”  category  and  consider  this  and  “tectonic”  as  our  only 
categories  of  interest. 

2.1.3  Availability  of  Regional  Phases 

The  REB  is  the  final  product  of  the  Prototype  IDC.  We  query  the  REB  database  to 
determine  the  number  of  regional  phases  that  are  available  for  each  subnetwork  defined  in 
Section  2.1.1.  We  use  the  following  constraints  to  limit  our  query  to  regional  events 
recorded  by  stations  in  the  subnetwork: 

•  Origin  time  between  January  1, 1995  and  June  1, 1996 

•  Associated  regional  phases  in  the  REB  (e.g.,  Pn.  Pg,  Sn,  Lg) 

•  Distance  less  than  20° 

•  Depth  less  than  40  km 

We  examine  the  quantity  of  regional  airivals,  the  availability  of  multiple-station 
associations,  and  the  tectonic  environment  for  each  subnetwork.  The  results  are 
summarized  in  Table  1.  Table  A1  in  Appendix  A  provides  more  detail  on  the  statistics  for 
each  station  in  each  subnetwork. 


Table  1:  SubNetworks  of  IDC  Primary  and  Auxiliary  Stations. 


SubNetwork 

Pn 

Pg 

Sn 

Lg 

Multiple- 

Station 

Events 

Tectonic 

Name 

Stations 

E.  Caucasus/Hindu 
Kush/ E,  Med.  Sea 

262 

27 

81 

35 

14 

Tectonic 

ABKT,  KBZ.  KVAR.  NIL,  BGIO 

Pacific  Rim 

1299 

60 

402 

93 

240 

Tectonic 

MJAR.OGS,  TSK 

Spain 

273 

23 

62 

67 

60 

Tectonic 

ESDC,  PAB 

Thailand 

265 

86 

23 

84 

0 

Tectonic 

CMAR 

Western  USA  / 
Western  Canada 

1150 

334 

227 

502 

267 

Tectonic 

PDAR.  PFO.  WALA,  BBB,  DUG,  EDM.  ELK, 

MNV,  NEW,  PGC,  PMB,  PNT,  RSSD,  TUC 

Alaska/NW  Canada 

1315 

154 

440 

438 

285 

Mixed 

WHY.  YKA.  NPO.  MBC,  DAWY,  DLBC,  INK 

Central  Africa 

53 

1 

29 

12 

0 

Mixed 

BGCA.  DBIC.  AAE 

South  America 

642 

12 

105 

70 

175 

Mixed 

BDFB.  CPUP,  LPAZ,  PLCA,  ROSC 

Antarctica 

39 

0 

22 

9 

0 

Stable 

MAW.  VNDA,  CSY 

Australia 

3683 

75 

2064 

339 

1028 

Stable 

ASAR.  STKA.  WOOL.  WRA.  ARMA.  CTA,  FITZ, 
FORT.  MEEK.  NWAO.  QIS,  RMQ,  TOO.  WARB, 
YOU 

Central  Asia 

197 

21 

71 

84 

19 

Stable 

ALFM,  BIT,  HIA,  NRI,  PDY,  ZAL,  ARU,  ULN, 

USK 

North  American 

Shield 

452 

34 

159 

159 

54 

Stable 

LBNH,  SCHQ,  ULM,  DRLN,  EEO,  EYMN,  FCC, 
FRB,  GAC,  LMN.  LMQ,  RES,  SADO,  TBO 

South  Africa 

48 

15 

31 

38 

9 

Stable 

BOSA,  1.SZ,  TSUM 

Southern  USA 

551 

60 

75 

190 

133 

Stable 

MIAR,  TXAR,  AI.Q,  TKL,  TUL 

Several  subnetworks  have  a  sufficient  number  of  regional  arrivals  for  our  attenuation  study 
(Australia,  Pacific  Rim,  Alaska/NW  Canada,  Western  USAAVestern  Canada,  South 
America,  Southern  USA,  and  North  America  Shield).  The  Central  Asia,  Thailand,  E. 
Caucasus/Hindu  Kush/E.  Med.  Sea,  and  Spain  subnetworks  may  have  enough  regional 
phases,  but  may  require  more  depending  on  the  geographic  distribution.  Other 
subnetworks  such  as  Central  Africa,  South  Africa,  and  Antarctica  do  not  appear  to  have 
enough  regional  phases  in  the  REB.  There  may  be  additional  regional  phases  recorded  at 
these  stations,  but  the  events  do  not  appear  in  the  REB  because  they  do  not  satisfy  the 
GSETT-3  event  confirmation  criteria.  Many  National  Data  Centers  (NDCs)  provide  event 
bulletins  based  on  local  networks  (called  Gamma  bulletins)  to  the  IDC.  As  described  in  the 
next  section,  these  local  events  can  be  used  to  supplement  our  data  sets. 

Table  1  shows  that  excellent  coverage  is  available  for  most  “stable”  regions,  and  several  of 
the  tectonic  regions.  Several  of  the  subnetworks  have  stations  in  both  stable  and  tectonic 
environments,  and  we  classify  these  subnetworks  as  “mixed.”  We  will  develop  attenuation 
models  for  each  subnetwork,  and  then  combine  data  from  several  subnetworks  to  derive 
generalized  attenuation  models  for  each  tectonic  environment  (e.g.,  tectonic  and  stable). 

We  may  change  the  subnetworks  as  we  accumulate  more  data  and  ground-truth 
information,  but  we  believe  that  the  14  subnetworks  in  Table  1  provide  the  best  possible 
data  sets  for  the  attenuation  study.  We  will  chose  our  final  10  subnetworks  from  this  list. 


2.2  Event  Selection 

We  search  event  bulletins  for  event  and  arrival  data  that  will  be  used  for  calculating  the 
attenuation  models.  We  use  the  same  selection  criteria  as  our  station  selection  criteria. 
Waveform  data  that  corresponds  to  our  selected  events  will  come  from  IDC  primary  and 
auxiliary  stations. 

We  also  compile  data  sets  that  will  be  used  to  evaluate  the  effectiveness  of  the  Lg/P 
discriminant.  These  are  likely  to  include  the  events  used  to  develop  the  attenuation  models, 
but  they  must  also  include  others.  The  attenuation  models  are  based  on  the  highest  quality 
data  that  are  available,  but  the  Lg/P  discriminant  must  be  tested  on  typical  events  in  the 
REB  (i.e.,  not  just  the  high-quality  events). 

2.2.1  Attenuation  Models 

The  REB  currently  includes  over  31,000  events  and  will  be  the  primary  source  of  data  for 
the  attenuation  models.  However,  this  is  supplemented  with  the  other  bulletins  whenever 
possible  to  fill  in  regions  with  sparse  coverage.  Twenty-one  National  Data  Centers  (NDCs) 
from  around  the  world  provide  event  bulletins  to  the  IDC.  The  compilation  of  these 
individual  NDC  bulletins  is  called  the  Gamma  bulletin,  and  we  use  this  bulletin  as  our 
main  supplemental  bulletin.  The  Gamma  bulletin  includes  many  local  events  that  are  not  in 
the  REB  due  to  strict  GSETT-3  event  confirmation  criteria,  the  small  size  of  local  events, 
and  sparse  local  IDC  network  coverage.  However,  some  detections  from  these  smaller 
events  may  have  been  recorded  by  IDC  primary  or  auxiliary  stations.  We  may  consider 
using  other  event  bulletins  such  as  the  Preliminary  Determination  of  Epicenters  (PDE) 
bulletin  produced  by  the  United  States  Geological  Survey  (USGS)  as  needed. 
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For  the  attenuation  models,  we  must  asses  the  ray-path  distribution  for  each  subnetwork. 
For  example,  Figure  4  shows  Pn,  Pg,  Sn,  and  Lg  ray-paths  for  our  Australian  subnetwork. 
We  use  these  maps  to  evaluate  whether  or  not  we  need  to  add  or  subtract  data  to  obtain 
adequate  coverage  for  a  region.  In  this  example,  there  are  many  more  Pn  phases  from 
events  originating  in  the  Java  Trench  to  the  north  than  from  events  to  the  south  or  east.  We 
may  want  to  remove  some  or  all  of  the  Java  Trench  events  to  avoid  biasing  our  results  for 
paths  from  that  source  region.  One  approach  would  be  to  only  use  data  from  Java  Trench 
events  with  high  signal-to-noise  ratio.  Alternately,  we  may  omit  events  from  the  Java 
Trench  so  that  our  attenuation  estimates  are  made  solely  from  events  within  the  continent 
of  Australia.  We  will  evaluate  which  strategy  is  most  effective. 

2. 2. 1.1  Supplemental  Bulletins 

Adding  events  and  arrivals  to  our  data  set  from  supplemental  bulletins  may  significantly 
improve  coverage  in  many  regions.  Table  2  lists  the  participating  NDCs  and  the  current 
number  of  events  in  the  Gamma  bulletin  that  were  reported  during  GSETT-3.  Figure  5 
shows  the  global  distribution  of  Gamma  events  that  might  be  used  in  our  inversion. 
However,  in  order  to  determine  whether  any  of  these  events  can  be  used  in  our  attenuation 
study,  we  must  develop  procedures  and  criteria  to  search  for  the  IDC  arrivals  that  are 
associated  with  events  in  the  Gamma  or  other  bulletins. 


Table  2:  Supplemental  NDC  Bulletins. 


Country 

National  Data  Center 
(NDC) 

Number  of  Events 

Australia 

AUS 

160 

Bulgaria 

BGR 

89 

Canada 

CAN 

475 

Switzerland 

CHE 

26 

China 

CHN 

3138 

Germany 

DEU 

147 

Denmark 

DMK 

3 

Spain 

ESP 

1059 

Finland 

HN 

1558 

France 

FRA 

281 

Great  Britain 

GBR 

34 

Hungary 

HUN 

13 

Israel 

ISR 

37 

Italy 

ITA 

488 

Japan 

JPN 

10110 

New  Zealand 

NZL 

3599 

Poland 

POL 

461 

Romania 

ROM 

255 

South  Africa  ■ 

RSA 

790 

Russia 

RUS 

1916 

United  States 

USA 

3007 
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Figure  4.  Example  from  the  Australian  subnetwork  showing  regional  events  recorded 
by  IDC  primary  and  auxiliary  stations  as  reported  in  the  Reviewed  Event  Bulletin. 
This  data  set  covers  the  time  period  of  January  1,  1995  to  June  1,  1996.  Many  more 
regional  phases  originate  from  events  to  the  north  than  from  events  to  the  south  or 
east,  (a)  Pn,  (b)  Sn^  (c)  Pg,  and  (d)  Lg  ray-paths. 
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Figure  4.  (Continued) 
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Figure  5.  Global  distribution  of  events  in  the  IDC  Gamma  bulletin. 


We  use  the  automated  association  software  at  the  IDC  to  associate  arrivals  in  the  REB  with 
events  in  the  Gamma  bulletin.  The  software  is  called  the  Global  Association  (GA)  system 
[LeBras  et  al,  1994].  It  uses  a  data-driven  exhaustive  grid-search  technique,  which  is 
based  on  the  generalized  beam  forming  method  introduced  by  Ringdal  and  Kvcema  [1989]. 

We  consider  arrivals  that  have  been  detected  by  the  automated  system  or  added  by  an 
analyst.  Using  arrivals  from  the  REB  in  a  5-minute  time  interval  around  the  origin  time  of 
Gamma  events,  GA  forms  event  hypotheses  over  a  regional  grid  around  the  subnetwork.  At 
each  grid  point,  GA  identifies  DRIVER  arrivals  and  possible  corroborating  regional 
arrivals.  A  DRIVER  is  an  arrival  at  one  of  a  limited  number  of  stations  in  the  network  that 
could  record  the  earliest  arrival  for  an  event  in  a  given  grid  cell  [LeBras  et  al.,  1994].  GA 
first  constrains  the  corroborating  arrivals  using  travel  time  and  slowness.  A  more  rigorous 
chi-square  statistical  test  is  applied  if  an  arrival  passes  this  initial  screening,  and  we  use  all 
available  features  of  the  arrival  (travel-time,  slowness,  azimuth,  and  amplitude).  GA  also 
allows  for  some  of  these  phases  to  be  identified  differently  than  the  initial  automated 
station  processing,  and  eliminates  outliers  based  on  a  chi-square  test  of  the  residuals  for 
each  arrival. 

After  these  tests  are  completed,  the  event  is  located  using  a  fixed  depth,  and  an  outlier 
analysis  is  applied.  Outliers  are  removed  if  necessary  and  the  location  is  refined.  GA  then 
applies  final  event  confirmation  criteria,  which  include  a  weighted-count  of  the  number  of 
defining  observations,  size  of  the  error  ellipse  major  axis,  and  the  probability  of  detection 
for  events  with  less  than  five  associations.  The  probability  of  detection  screens  event 
hypotheses  that  were  detected  by  an  unlikely  combination  of  stations. 

As  shown  in  Figure  4a,  the  Australian  region  is  dominated  by  events  to  the  north,  and 
additions  to  the  data  set  from  the  Gamma  bulletin  may  improve  our  coverage  within  the 
continent.  We  applied  GA  to  the  160  events  provided  by  the  Australian  NDC.  Of  these,  101 
are  not  in  the  REB.  For  these  non-REB  events,  we  found  44  arrivals  in  the  REB  that  appear 
to  be  associated  with  17  regional  events.  Figure  6  shows  the  raypaths  for  the  44  arrivals. 
Adding  these  phases  will  contribute  to  our  coverage  of  central  Australia.  The  remaining  84 
events  could  not  be  formed  from  IDC  arrivals  due  to  insufficient  data  from  our  primary  and 
auxiliary  stations. 

We  compared  the  analyst-verified  REB  associations  to  those  that  we  obtained  for  the 
events  in  common  with  both  the  Gamma  bulletin  and  the  REB.  When  comparing  the  total 
number  of  phases,  we  obtain  a  high  percentage  of  associated  phases  (93%)  and  defining 
phases  (80%)  that  are  in  the  REB.  Of  our  defining  phases,  82%  exactly  match  the  REB 
defining  phases.  Thus,  we  have  confidence  that  our  procedure  is  valid  since  it  produced 
very  similar  association  sets  and  locations  for  events  in  common  with  the  REB. 

2.2.2  Lg/P  Discriminant 

We  will  assemble  two  test  data  sets  to  evaluate  the  effectiveness  of  the  Lg/P  discriminant 
for  each  subnetwork.  The  first  data  set  will  consist  of  events  whose  identifications  are 
known  with  high  confidence,  and  we  call  it  the  ground-truth  data  set.  It  will  be  assembled 
by  searching  event  bulletins  that  identify  source  type.  The  Gamma  bulletin  will  be  used  for 
the  ground-truth  data  set,  since  many  events  are  identified  (mine  blast,  explosion. 
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Figure  6.  Ray  paths  for  arrivals  associated  with  the  Gamma  bulletin  provided  by  the 
Australian  National  Data  Center  (NDC).  (a)  Pn,  (b)  Sn,  (c)  Pg,  and  (d)  Lg. 
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Figure  6.  (Continued) 
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earthquake,  etc.).  However,  we  expect  that  we  will  need  to  supplement  this  with  other 
ground-truth  information  from  participating  NDCs. 

The  second  data  set  will  include  a  large  number  of  events  (>200  per  subnetwork)  intended 
to  represent  the  type  of  events  that  are  routinely  identified  by  a  seismic  monitoring  system. 
In  general,  we  will  not  have  ground-truth  information  for  these  events.  Instead,  we  will 
rely  on  knowledge  of  the  local  natural  and  industrial  seismicity  to  qualitatively  evaluate  the 
effectiveness  of  the  Lg/P  discriminant.  We  will  use  the  approach  developed  by  Sereno  and 
Wahl  [1993]  for  evaluating  the  Event  Identification  System  in  ARPA’s  Intelligent 
Monitoring  System  (IMS).  They  used  knowledge  of  the  local  natural  and  industrial 
seismicity  to  qualitatively  evaluate  the  effectiveness  of  discriminants.  Figure  7  shows  an 
example  of  results  using  this  approach.  All  of  the  events  that  were  identified  as  mine  blasts 
occurred  in  areas  of  active  mining.  Most  of  the  events  that  were  identified  as  earthquakes 
are  in  areas  with  the  highest  natural  seismicity  in  this  region  (the  southwest  coast  of 
Norway  and  the  Mid- Atlantic  Ridge).  Many  of  the  events  that  were  identified  as 
explosions  are  in  active  mining  areas,  but  the  case-based  approach  used  in  the  IMS  failed 
to  associate  them  with  a  known  mine.  Several  of  the  offshore  clusters  are  known  to  be 
explosions  from  military  exercises.  Some  of  the  events  are  obviously  identified  incorrectly 
by  the  automated  system  (such  as  the  events  that  were  identified  as  explosions  on  the  Mid- 
Atlantic  Ridge).  However,  this  example  shows  how  we  can  gain  valuable  information 
regarding  the  performance  of  a  discriminant  (or  set  of  discriminants)  using  the  general 
properties  of  a  large  sample  of  events  without  ground-truth  information. 
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Earthquakes 


Mine  Blasts  Explosions  (ML  >  2.0) 


Figure  7.  These  maps  show  the  location  of  events  that  were  identified  as 
earthquakes,  mine  blasts  and  explosions  by  the  IMS  Event  Identification  System. 
Ground-truth  information  is  not  available  for  these  events,  but  it  is  clear  that  most 
of  the  events  are  identified  correctly.  The  number  of  events  is  listed  in  the  lower 
right  comer  of  each  map. 
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3.  Attenuation  Models 

We  will  estimate  frequency-dependent  attenuation  of  regional  phases  recorded  by  up  to  10 
IDC  subnetworks.  We  measure  amplitudes  of  regional  phases,  and  develop  software  to 
invert  them  for  source  and  attenuation  models.  Figure  8  shows  the  major  software 
components.  DFX  was  developed  for  the  PIDC  to  perform  seismic  signal  detection  and 
feature  extraction.  It  will  be  used  in  this  study  to  compute  regional  wave  amplitudes  in  the 
time-  and  frequency-domains.  Ampinv  is  a  new  program  being  developed  under  this 
contract  to  invert  the  regional  wave  amplitudes  for  source  and  attenuation  parameters. 


Figure  8.  DFX  and  Ampinv  are  the  two  major  software  components  for  this 
project.  Both  programs  interface  with  the  IDC's  Data  Management  System. 


3.1  Amplitude  Measurements 

The  software  program  DFX  (Detection  and  Feature  extraction)  was  developed  for  the 
PIDC  to  perform  seismic  signal  detection  and  feature  extraction  [Wahl,  1996].  For  this 
project,  we  added  the  capability  for  DFX  to  compute  amplitude  spectra. 

We  investigate  three  types  of  time-domain  amplitude  measurements  for  both  3-component 
and  array  stations.  These  are  root  mean-square  {rms),  short-term  average  {stav),  and  the 
average  of  absolute  amplitude  in  the  time  window  after  removing  the  rms.  The  rms 
measurements  have  been  used  in  regional  magnitude  estimation  and  explosion  yield 
studies  because  they  provide  stable  amplitude  estimates  for  regional  phasesje.g.,  Ringdal, 
1983;  Nuttli,  1986fl,h;  Israelsson,  1992].  We  evaluate  the  different  amplitude  measures  by 
comparing  the  variance  with  respect  to  our  derived  attenuation  models.  Our  preferred 
amplitude  measurement  for  each  phase  is  the  one  with  the  lowest  variance. 

Measurements  are  made  on  both  coherent  and  incoherent  beams  in  up  to  seven  overlapping 
frequency  bands  (each  one  octave  wide)  from  0.5  Hz  to  16.0  Hz.  We  use  10-second  signal 
windows  for  Pn,  Pg,  and  Sn  starting  0.3  seconds  before  the  analyst-verified  onset  time  in 
the  REB.  A  group  velocity  window  of  3.0-3. 6  km/sec  is  used  for  Lg.  Noise  windows  are 
5.0  seconds  and  start  5.5  seconds  prior  to  the  signal  window.  We  calculate  spectra  by 
demeaning  and  cosine-tapering  the  windowed  data  and  applying  an  FFT.  For  array 
stations,  the  spectra  are  array-averaged  using  the  method  of  Bache  et  al.  [1985]. 
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Measurements  will  be  made  on  vertical  data  channels,  but  we  may  evaluate  other  channels 
in  the  future  (e.g,  radial,  transverse,  and  combined  horizontal). 

Frequencies  with  adequate  signal-to-noise  ratio  (snr)  are  determined  for  each  phase. 
Determination  of  the  frequency  limits  is  bounded  by  the  instrument  response  at  the 
recording  station  and  limits  provided  by  the  user.  The  snr  is  required  whether  time-  or 
frequency-domain  amplitude  measurements  are  used.  For  the  time-domain,  we  simply 
compare  the  snr  to  the  user-specified  threshold.  For  the  frequency-domain  (i.e.  amplitude 
spectra),  we  smooth  the  signal  and  noise  curves  and  resample  them  at  a  coarse  sampling 
rate  specified  by  a  user-parameter.  The  snr  of  each  sample  must  exceed  the  threshold 
before  that  sample  can  be  used  in  the  inversion. 

We  have  computed  approximately  6000  amplitude  measurements  for  regional  phases 
recorded  at  Australian  stations.  As  an  example.  Figure  9a  shows  a  regional  seismogram 
recorded  by  the  Alice  Springs  array  (ASAR)  in  Australia  from  an  Ml=3.6  event  in  south- 
central  Australia  at  a  distance  of  about  500  km.  The  analyst-reviewed  onset  times  are 
superimposed.  We  use  a  group  velocity  window  for  Lg  based  on  epicentral  distance  and 
estimates  of  the  event  origin  time.  Figure  9b  shows  that  this  time  window  captures  most  of 
the  analyst-reviewed  Lg  onset  times.  It  shows  arrival  times  for  210  Lg  arrivals  recorded  by 
stations  in  Australia.  Amplitude  spectra  are  computed  using  the  same  time  windows  as 
those  used  for  time-domain  measurements.  Figure  9c  shows  smoothed  and  instrument- 
corrected  signal  and  noise  spectra  for  the  Pn  arrival  in  Figure  9a.  We  use  these  spectra  to 
determine  which  frequencies  to  include  in  the  inversion. 

3.2  Inversion  Method  and  Software 

We  use  the  method  of  estimating  attenuation  developed  by  Sereno  [1990].  Briefly,  the 
frequency-dependent  (either  time-domain  or  spectral)  amplitude  of  the  kth  wave  recorded 
at  the  ith  station  from  the  jth  source  is  parameterized  as: 

if)  =  logA"  if)  +  S,  ( A„,/)  +  (1) 

where  is  the  amplitude  at  a  reference  distance  Aq,  Aq,  f)  is  the  attenuation 

from  the  reference  distance  to  the  epicentral  distance  Ay,  and  is  a  station  correction.  The 
amplitude  at  the  reference  distance  is  expressed  in  terms  of  the  material  properties  at  the 
source  and  the  receiver,  source  parameters  such  as  the  seismic  moment,  the  shape  of  the 
source  spectrum,  and  a  wave-dependent  excitation  factor.  We  use  the  Mueller  and  Murphy 
[1971]  and  Brune  [1970,  1971]  models  to  parameterize  the  shape  of  the  source  function. 

The  attenuation  is  parameterized  in  terms  of  a  power-law  distance  dependence  with  a 
frequency-dependent  exponent: 

^k  ^ 

where  the  first  term  accounts  for  attenuation  from  the  source  to  the  reference  distance,  and 
the  second  term  is  the  total  attenuation  from  the  reference  distance  to  the  epicentral 
distance.  The  total  attenuation  includes  geometrical  spreading,  scattering  and  anelasticity. 
It  is  difficult  to  separate  these  terms  since  the  geometrical  spreading  of  regional  phases  is  a 
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(a) 


Si 


(b)  Lg  Group  Velocity  Window 


Figure  9.  (a)  Example  of  a  waveform  from  one  of  the  elements  of  the  ASAR  array.  The 
majority  of  the  Lg  energy  is  included  in  the  3.0-3.6  km/s  group  velocity  window,  (jb) 
Arrival  times  for  210  Lg  phases  in  the  IDC  database  recorded  by  stations  in  Australia. 
Most  arrivals  are  included  in  the  3.0-3.6  km/ s  group  velocity  window,  (c)  Signal  and 
noise  amplitude  spectra  for  the  Pn  arrival  shown  in  (a). 
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Figure  9.  (Continued) 

complicated  function  of  the  crustal  and  upper  mantle  velocity  structure  [e.g.,  Sereno  and 
Given,  1990],  Fortunately,  it  is  not  important  to  separate  these  terms  for  application  to  the 
Lg/P  discriminant.  The  exponent,  njjif),  is  parameterized  as  a  linear  function  of  frequency. 

The  system  of  equations  governing  the  relationship  between  the  data  (log  amplitudes)  and 
model  parameters  is  formulated  by  subtracting  theoretical  data  computed  from  an  assumed 
starting  model  from  the  observed  data,  and  solving  iteratively  for  the  model  perturbations 
that  minimize  the  data  residuals  in  a  least-squares  sense.  The  inversion  process  is  described 
in  detail  in  Figure  10.  We  calculate  the  data  residual  {dres)  between  the  observed  data 
{dobs)  and  data  predicted  {dpre)  based  on  the  starting  model  (M).  Next,  we  create  the 
partial  derivative  matrix  containing  the  system  of  equations  (G)  describing  our  problem 
and  then  decompose  it  by  a  singular  value  decomposition  (SVD)  as  outlined  by  Stewart 
and  Leyk  [1994].  The  decomposition  yields  eigenvector  matrices  for  the  model  (V)  and 
data  (U),  and  also  calculates  a  matrix  of  singular  values  (A).  Damping  is  applied  as  a 
function  of  the  condition  number  of  the  singular  values  as  required.  The  degree  of  damping 
is  controlled  by  the  magnitude  of  the  condition  number.  The  model  perturbation,  5m,  is 
determined  as: 

5m  =  V  *  *  dres  (3) 

where  V  is  the  matrix  of  model  eigenvectors,  A“^  is  the  inverse  of  singular  values  matrix, 
is  the  transpose  of  the  matrix  of  data  eigenvectors,  and  dres  is  the  data  residuals.  We 


21 


add  the  model  perturbation  to  the  model  and  proceed  with  the  next  iteration.  This  continues 
until  convergence.  Our  convergence  criteria  are  based  on  the  maximum  number  of 
iterations,  the  condition  number  of  the  singular  values  matrix  (A),  and  a  data  variance 
threshold. 


Start 


Figure  10.  The  inversion  process  receives  an  observed  data  set  (dobs)  and  a 
starting-model  (M).  A  singular  value  decomposition  (SVD)  method  is  used  to 
obtain  data  eigenvector  (U),  model  eigenvector  (V),  and  singular  value  (A) 
matrices.  Adjustments  are  made  to  the  model  (5m)  until  the  convergence  criteria 
are  met.  Then,  the  final  model  (M'),  and  model  covariance  (Cov),  model 
resolution  (R),  and  data  resolution  (N)  matrices  are  returned. 


Ampinv  is  a  new  software  program  that  we  are  developing  to  implement  the  inversion 
method  described  above  (Figure  11).  It  will  invert  frequency-dependent  amplitudes  for 
source  and  attenuation  parameters  such  as  the  long-period  source  level,  a  scaling  parameter 
relating  comer  frequency  to  long-period  source  level,  phase-dependent  excitation  factors, 
anelastic  attenuation  from  the  source  to  a  referenee  distance,  total  attenuation  from  the 
reference  distance  to  receiving  station,  and  a  phase-dependent  station  correction  [Sereno, 
1990].  Ampinv  also  generates  synthetic  amplitudes  with  or  without  random  noise.  We  will 
use  synthetic  amplitudes  to  validate  Ampinv,  experiment  with  user  parameters  (e.g., 
convergence  criteria),  and  evaluate  the  resolution  of  source  and  attenuation  parameters  for 
each  data  set.  Ampinv  obtains  event  and  amplitude  information  directly  from  the  DMS 
using  a  conneetion  through  our  Generic  Database  Interface  [Anderson  et  al,  1994].  The 
time-domain  amplitudes  are  stored  in  an  Oracle™  Relational  Database  Management 
System  [Swanger  et  al,  1993],  and  the  amplitude  spectra  are  stored  as  Unix  files  on  disk. 
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Figure  11.  Amplnv  process  flow  diagram.  Thin  lines  link  main  functional  imits 
within  Amplnv.  Thick  lines  show  interaction  with  user  parameters,  the  Data 
Management  System  (DMS),  and  output  model  file. 


The  User’s  Guide  for  Amplnv  gives  a  detailed  description  of  the  user  parameters, 
algorithms,  and  data  stractures  [Jenkins,  1996].  User  parameters  including  database, 
signal-to-noise  limits,  phase  dependent  thresholds,  starting  model  and  other  parameters 
that  control  Amplnv  are  normally  stored  in  a  parameter  (or  par)  file,  but  can  also  be 
specified  on  the  command-line.  Amplnv  reads  user  input,  queries  the  IDC  Data 
Management  System  (DMS)  for  input  data,  computes  synthetic  amplitude  if  desired, 
determines  usable  frequency  bands  for  each  arrival,  inverts  for  source  and  attenuation 
parameters  using  the  inversion  process  in  Figure  10,  and  writes  an  output  model. 
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4.  Summary  and  Future  Work 

During  the  first  year  of  this  two-year  project,  we  developed  all  software  needed  for  this 
study,  investigated  alternative  amplitude  measures,  and  compiled  a  data  set  for  an 
Australian  subnetwork.  During  the  second  year,  we  will  compile  the  data  sets  for  the  other 
IDC  subnetworks,  estimate  regional  wave  attenuation  using  Ampinv,  and  evaluate  the 
accuracy  and  transportability  of  the  Lg/P  discriminant.  We  will  also  obtain  “ground-truth” 
identification  for  as  many  of  the  events  as  possible,  and  we  will  use  knowledge  of  the  local 
natural  and  industrial  seismicity  when  this  information  is  not  available. 
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Appendix  A:  IDC  Primary  Station  Areas 

Table  A1  lists  the  network  name,  station  type  (e.g.,  Pri=primary,  Aux=auxiliary),  station 
name,  number  of  regional  phases,  number  of  regional  and  multiple-station  events,  and  the 
tectonic  environment  of  each  member  of  the  network. 


Table  Al:  IDC  Primary  Station  Areas 


Subnetwork 

Sta 

Type 

Station  Name  /  IDC  Name 

Pn 

Pg 

Sn 

Lg 

Num. 

Events 

Multiple- 

Station 

Events 

Tectonic 

Name 

Alaska/NW 

Canada 

Pri 

Whitehorse,  CanadaAVHY 

305 

39 

104 

101 

309 

285 

Paleo/Mesozoic 

Yellowknife,  CanadaA’KA 

279 

13 

78 

285 

205 

Archean 

North  Pole,  Alaska/NPO 

315 

38 

103 

101 

319 

243 

Cenozoic 

Mould  Bay,  CanadaAlBC 

135 

4 

56 

14 

135 

114 

Paleo/Mesozoic 

Aux 

Dawson  City,  Canada/DAWY 

190 

37 

54 

72 

195 

188 

Paleo/Mesozoic 

Dease  Lake,  Canada/DLBC 

108 

12 

33 

42 

109 

104 

Paleo/Mesozoic 

Inuvik.  Canada/INK 

118 

15 

58 

44 

118 

118 

Paleo/Mesozoic 

Total 

1450 

158 

496 

452 

Antarctica 

Pri 

Mawson,  Antarctica/MAW 

19 

0 

17 

0 

19 

0 

Archean 

Vanda,  Antarctica/VNDA 

13 

0 

5 

2 

13 

0 

Paleo/Mesozoic 

Aux 

Casey  Station,  Antarctica/CSY 

7 

0 

0 

7 

7 

0 

Archean 

Total 

39 

0 

22 

9 

Australia 

Pri 

Alice  Springs,  Au.stralia/ASAR 

812 

15 

642 

97 

823 

819 

Proterozoic 

Stephens  Creek,  Aus/STKA 

106 

20 

64 

53 

106 

102 

Proterozoic 

Woolibar,  AustraliaAVOOL 

78 

5 

41 

18 

80 

75 

Archean 

Warramunga,  AustraliaAVRA 

1560 

21 

961 

97 

1581 

1028 

Proterozoic 

Aux 

Armidale,  Australia/ARMA 

38 

0 

2 

1 

38 

13 

Paleo/Mesozoic 

Charters  Towers,  Australia/CTA 

184 

1 

22 

13 

188 

78 

Paleo/Mesozoic 

Fitzroy  Crossing,  Aus/FITZ 

5 

0 

2 

0 

5 

4 

Archean 

Forrest,  Australia/FORT 

21 

0 

9 

6 

23 

22 

Archean 

Meekathara,  Australia/MEEK 

126 

1 

63 

4 

127 

99 

Archean 

Narrogin,  Australia/NWAO 

II 

2 

4 

5 

14 

11 

Archean 

Mt  Isa,  Australia/QIS 

455 

0 

138 

15 

457 

388 

Proterozoic 

Roma,  Australia/RMQ 

67 

8 

0 

6 

66 

35 

Toolangi,  Australia/TOO 

18 

0 

7 

4 

18 

11 

Paleo/Mesozoic 

Warburton,  Australia/WARB 

197 

1 

109 

3 

198 

194 

Proterozoic 

Young,  Australia/YOU 

5 

1 

3 

0 

6 

5 

Paleo/Mesozoic 

Total 

3683 

75 

2064 

Table  Al:  IDC  Primary  Station  Areas 


Subnetwork 

Sta 

Type 

Station  Name  /  IDC  Name 

Pn 

Pg 

Sn 

Lg 

Num. 

Events 

Multiple- 

Station 

Events 

Tectonic 

Name 

Central  Africa 

Pri 

Bogoin,  C.  African  Rep/BGCA 

18 

0 

16 

11 

18 

0 

Archean 

Dimbroko,  Ivory  Coast/DBIC 

30 

■1 

13 

1 

30 

0 

Archean 

Aux 

Addis  Ababa,  Ethiopia/AAE 

5 

0 

0 

0 

5 

0 

Paleo/Mesozoic 

Total 

53 

IHi 

29 

12 

Central  Asia 

Pri 

Mangolia/ALFM 

0 

0 

0 

0 

0 

0 

Paleo/Mesozoic 

Baijiatuan,  China/BJT 

24 

3 

4 

6 

24 

6 

Archean 

Hailar,  China/HIA 

3 

0 

0 

0 

3 

1 

Paleo/Mesozoic 

Norilsk,  Russia/NRI 

14 

0 

5 

4 

15 

14 

Archean 

Peleduy,  Russia/PDY 

62 

6 

23 

33 

64 

19 

Archean 

Zalesovo,  Russia/ZAL 

59 

■B 

20 

23 

61 

17 

Paleo/Mesozoic 

Aux 

Arti,  Russia/ARU 

23 

■1 

14 

9 

26 

6 

Paleo/Mesozoic 

Ussuriisk,  Russia/USK 

0 

0 

0 

0 

0 

0 

Paleo/Mesozoic 

Ulaanbaatar,  Mongolia/ULN 

12 

4 

5 

9 

12 

6 

Paleo/Mesozoic 

Total 

197 

21 

71 

84 

East  Caucasus/ 
Hindu  Kush/ 

E.ast  Med.  Sea 

Pri 

Alibek,  Turkmenistan/ABKT 

32 

2 

13 

■ 

32 

10 

Cenozoic 

Khabaz,  Russia/KBZ 

96 

8 

12 

10 

97 

24 

Tectonic 

Aux 

Kislovodsk  Array,  Russia/KVAR 

24 

3 

5 

2 

25 

10 

Tectonic 

Bar  Giyora,  Israel/BGIO 

53 

3 

21 

6 

55 

11 

Proterozoic 

Nilore,  Pakistan/NIL 

57 

11 

30 

13 

59 

4 

Proterozoic 

Total 

262 

27 

81 

35 

North 

American  Shield 

Pri 

Lisbon,  New  Hampshire/LBNH 

15 

3 

6 

15 

18 

18 

Proterozoic 

Schefferville,  Canada/SCHQ 

31 

0 

12 

9 

31 

17 

Archean 

Lac  du  Bonnet,  Canada/ULM 

121 

13 

17 

54 

126 

48 

Archean 

Aux 

Deer  Lake,  Canada/DRLN 

23 

0 

6 

0 

23 

4 

Paleo/Mesozoic 

Eldee,  Canada/EEO 

3 

1 

1 

3 

3 

3 

Archean 

Ely,  Minnesota/EYMN 

6 

2 

4 

5 

9 

7 

Archean 

Fort  Churchill,  Canada/FCC 

56 

1 

20 

26 

59 

40 

Archean 

Iqaluit,  Canada/FRB 

8 

0 

6 

4 

8 

8 

Archean 

Glen  Almond,  Canada/GAC 

■01 

3 

4 

5 

6 

6 

Proterozoic 

Caledonia  Mtn.,  Canada/LMN 

2 

0 

1 

1 

2 

2 

Paleo/Mesozoic 

La  Malbaie,  Canada/LMQ 

5 

0 

5 

5 

8 

8 

Proterozoic 

29 


Table  Al:  IDC  Primary  Station  Areas 


Subnetwork 


Sta 

Type 

Station  Name  /  IDC  Name 

Pn 

Aux 

Resolute  Bay,  Canada/RES 

28 

Sadowa,  Canada/SADO 

13 

Thunder  Bay,  Canada/TBO 

2 

Num.  Tectonic 

Events  Name 

Events 


317  30  103  145 


Pacific  Rim 

Pri 

Matsushiro  Array,  Japan/MJAR 

956 

46 

244 

71 

Aux 

Ogasawara,  Japan/OGS 

165 

6 

53 

7 

Tsukuba,  Japan/TSK 

178 

8 

105 

15 

South  Africa  Pri 
Aux 


South  America  Pri 


Boshof,  South  Africa/BOSA 


Lusaka,  Zambia/LSZ 


Tsumeb,  Namibia/TSUM 


12  105 


Southern  USA  Pri 


75  190 


Thailand  Pri 


Total 


Chiang  Mai,  Thailand/CMAR 


Total 


Paleo/Mesozoic 


1 3  Proterozoic 


3  Archean 


240  Tectonic 


133  Tectonic 


177  Old  Ocean 


Pri 

Brasilia,  Brazil/BDFB 

4 

0 

0 

0 

4 

3 

Paleo/Mesozoic 

Villa  Florida,  Paraguay/CPUP 

188 

0 

33 

20 

188 

164 

Proterozoic 

La  Paz,  Bolivia/LPAZ 

218 

3 

34 

15 

217 

73 

Paleo/Mesozoic 

Paso  Flores,  Argentina/PLCA 

232 

9 

38 

35 

234 

175 

Tectonic 

El  Rosal,  Colombia/ROSC 

0 

0 

0 

0 

0 

0 

Tectonic 

Pri 

Mount  Ida,  Arkansas/MIAR 

108 

8 

22 

33 

f  109 

89 

Proterozoic 

TXAR  Array,  Texas/TXAR 

321 

25 

31 

117 

323 

133 

Paleo/Mesozoic 

Aux 

Albuquerque,  N  M/ALQ 

71 

18 

10 

23 

74 

60 

Proterozoic 

Tuckaleechee  Caverns,  TN/TKL 

9 

6 

4 

9 

13 

9 

Proterozoic 

Tulsa,  Oklahoma/TUL 

42 

3 

8 

8 

46 

44 

Proterozoic 

Spain 

Pri 

Sonseca  Array,  Spain/ESDC 

211 

16 

41 

47 

208 

60 

Paleo/Mesozoic 

Aux 

San  Pablo,  Spain/PAB 

62 

7 

21 

20 

62 

60 

Paleo/Mesozoic 

Western  USA/ 
Western 

Canada 

Pri 

Pinedale  Array,  WY/PDAR 

292 

77 

38 

109 

305 

267 

Paleo/Mesozoic 

Pinon  Flat,  Califomia/PFO 

121 

49 

24 

52 

137 

128 

Tectonic 

Table  Al:  IDC  Primary  Station  Areas 


Subnetwork 

Sta 

Type 

Station  Name  /  IDC  Name 

Pn 

Pg 

Sn 

Lg 

Num. 

Events 

Multiple- 

Station 

Events 

Tectonic 

Name 

Pri 

Waterton  Lakes,  CanadaAVALA 

202 

52 

44 

108 

233 

207 

Paleo/Mesozoic 

Aux 

Bella  Bella,  Canada/BBB 

42 

8 

24 

13 

47 

37 

Paleo/Mesozoic 

Dugway,  Utah/DUG 

64 

25 

6 

31 

71 

71 

Tectonic 

Edmonton,  Canada/EDM 

66 

10 

10 

23 

72 

71 

Paleo/Mesozoic 

Elko,  Nevada/ELK 

25 

10 

1 

13 

25 

24 

Paleo/Mesozoic 

Mina,  Nevada/MNV 

46 

30 

12 

26 

55 

55 

Paleo/Mesozoic 

Newport,  Washington/NEW 

35 

9 

7 

19 

38 

38 

Paleo/Mesozoic 

Pacific  Geo.  Ctr.,  Canada/PGC 

28 

IB 

13 

6 

31 

31 

Paleo/Mesozoic 

Pemberton,  Canada/PMB 

65 

10 

17 

23 

66 

61 

Paleo/Mesozoic 

Penticton,  Canada/PNT 

65 

12 

17 

24 

66 

66 

Paleo/Mesozoic 

Black  Hills,  SD/RSSD 

21 

20 

3 

22 

40 

38 

Archean 

Tucson,  Arizona/TUC 

78 

18 

11 

33 

86 

72 

Paleo/Mesozoic 

Total 

1150 

334 

227 

502 

31 
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